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Abstract 

/  - 


y  Although  many  of  the  rules  for  detecting  end  dealing  with 
initialization  bias  in  computer  simulation  experiments  are  easy  to 

-r?. 

understand  and  implement,  they  are  nonetheless  heuristic.  -The  current- 

paper  ur.es  the  theory  of  standardized  time  series  to  construct  tests 

'  9"" 

which  (under  certain  conditions)  detect  "signif icant"7  initialization 
bias  in  a  process.  Previous  tests  for  initialization  bias  can  be 
viewed  as  special  cases  of  the  general  family  of  tests  to  be  presented 
here. 
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INTRODUCTION 


When  a  system  is  simulated,  it  is  sometimes  difficult  to  select 
appropriate  initial  conditions  to  drive  the  simulation.  Me  often  wish 
to  initialize  the  process  with  "typical"  system  values;  however,  if 
there  is  little  a  priori  knowledge  of  the  process  available,  we  might 
initialize  the  system  in  a  state  which  has  a  very  low  probability  of 
occurring.  This  initialization  bias  (or  initial  transient)  can  then 
result  in  incorrect  conclusions  on  the  part  of  the  experimenter. 

A  common  way  of  dealing  with  the  initial  transient  is  to  take  a 
very  large  (possibly  wasteful)  number  of  observations  -  large  enough 
so  that  the  initialization  effects  are  overwhelmed.  Perhaps  a  better 
method  of  counteracting  the  bias  problem  is  simply  to  delete 
(truncate)  a  portion  of  the  output  from  the  beginning  of  the 
simulation  run.  The  experimenter  would  then  hope  that  the  offending 
biased  observations  had  been  eliminated.  Unfortunately,  if  the  output 
is  truncated  too  early,  then  significant  initialization  bias  might 
still  be  present.  If  it  is  truncated  too  late,  then  "good" 
observations  are  lost.  Ccf.  Snell  and  Schruben  (1979). 1 

Although  many  of  the  rules  for  detecting  and  dealing  with 
initialization  bias  are  easy  to  understand  and  implement,  they  are 
nonetheless  heuristic  [see,  e.g.,  the  surveys  by  Wilson  and  Pritsker 
(1978)  and  Schruben  and  Goldsman  (1984)1.  The  current  paper  uses  the 
theory  of  standardized  time  series  to  construct  tests  which  (under 
certain  conditions)  detect  "significant"  initialization  bias  in  a 
process.  Previous  tests  from  Schruben  (1982)  and  Schruben,  Singh,  and 
Tierney  (1983)  can  be  viewed  as  special  cases  of  the  general  family  of 
tests  to  be  presented  here. 

This  paper  is  organized  as  follows.  Background  material  is 
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provided  in  Section  2.  Our  new  teste  ere  introduced  in  Section  3. 
Section  4  is  concerned  with  power  calculations  f or  the  new  tests. 


2.  PRELIMINARIES 

2.1  Some  Standardized  Time  Series  Results 

Consider  the  stochastic  process  X ,,...,  X  .  For  jBl,...,m,  let 

l  (it 

X.  e  ,  X.  and  S.eX  -  X.;  also,  define  e  lim  mVar(X  ). 

j  j  l  j  m  J*  m-»®  m 

I mt I S  *  .  . 

The  standardized  time  series  is  T  (t)  s  - ; - -  ,  0  <  t  <  1,  where 

m  eVm  *  —  — 

l.|  is  the  greatest  integer  function. 


Suppose  that  X^,...,Xm  is  a  sequence  of  stationary,  finite 
variance  random  variables.  In  Schruben  (1983),  it  is  shown  that  under 
rather  mild  assumptions,  "^(t)  5  B^  as  m  -»  »,  where  B^  is  a  standard 
Brownian  bridge  process.  Further,  Tm(t)  is  asymptotically  independent 

of  mX  . 

m 

Let  us  divide  the  stationary,  finite  variance  series  X^,...,Xn 
into  b  adjacent  batches ,  each  consisting  of  m  X^'s  (n=bm) ;  the  random 

variables  X(i_1)m+1 *X(i_1)m+2 . X-m  comprise  batch  i,  i=l,...,b.  If 

m  is  large  enough,  we  can  treat  the  batches  as  if  they  were 
(approximately)  independent.  Each  individual  batch  can  then  be 
standardized;  this  yields  b  standardized  time  series  which  are 
approximately  independent  Brownian  bridges. 

For  i=l,...,b  and  j— l,...,m,  define  the  random  variables: 


r.  .  t  i  yJ  ,  X,. 

1,J  J  wp»l  (l-l)l 


(cumulative  averages). 


(Observe  that  X.  is  the  random  variable  corresponding  to  the  i-th 

l  ,m 


batched  mean . ) 


V  y 


y 


X  e  3T  "  —  I”  «  (grand  mean)  , 

n  b  t>i-l  i  ,m  n  «-i*l  1  "  ? 


1  cn 


S.  *  -  X. 

i  ,  j  i  ,m  i«J 


K.  e  argmax .  CkS.  .>, 

X  K  1  f  K 


S.  e  K.S.  £  .  and 
i  ii ,K. ’ 

l 


a.  s  ym  js.  .. 

1  fc«J«l  1,J 


Theorem:  We  have  the  -fallowing  collection  o-f  estimators  -for  a  t 

(0)  Classical  batched  means  estimator: 


..  0,b 

vo,b  £  nr  >  "here 


°n  h  1  *  ?  i  P  -  i  Jb  .  J  .  )2  2  »2*2(b-l>,  b  >  1 

0,b  ti*!1-  ifm  b  fcj8*!  j,mJ  * 


< 1 )  Area  estimator: 

vl,b  s  -E^  ’  "hBre 


12  cb  ~2  D  2  2..  ,  .  N  . 
),  .  A.  •*  c  X  lb)  ,  b  >  1. 
faxal  i  — 


lib  ~  3 

m  — m 


(2)  Combined  classical -area  estimator: 
Q 

v2,b  5  SE=T  ’  "hBrB 


Q0  .  s  Q.  .  +  Q,  .  5  <72X2  (2b-l )  ,  b  >  1, 

V|  D  1  y  b 


'  W  •  V> L  -X> 


i  ’r  .'C  C  i  C-.. 
/  r<  •  r  ■  >  .  k 

v  '■f  l»  i 


The  above  variance  estimators  Mill  be  used  in  this  paper  to  test 
output  -from  a  stochastic  process  Tor  initialization  bias. 

Remarks  Other  variance  estimators  arise  from  spectral  methods 
CHeidelberger  and  Welch  (1981)1  and  ARMA  time  series  modelling  methods 
C Fishman  (1971,1973,1978)  and  Andrews  and  Schriber  (1982)1;  the 
current  paper  will  not  concentrate  specifically  on  these  estimators. 


2.2  Previous  Initialization  Bias  Tests 
2.2.1  Motivation 

Suppose  that  we  model  X^,...,X  as  ,  i- l,...,n,  where 

ECXil  =  for  all  i,  and  the  's  are  stationary.  We  say  that  no 
initialization  bias  is  present  if  =  u,  say,  for  all  i.  Otherwise, 
bias  is  present.  Initialization  bias  might  exist  for  higher  order 
moments,  but  this  case  is  not  considered  Csee  Schruben  (1981)3. 

(ui>  is  the  transient  mean  process.  Figure  1(a)  is 
representative  of  so-called  negative  initialization  bias  (as  might  be 


u 

encountered  in  a  queue-length  process  when  starting  a  system  empty  and 
idle).  Figure  1 (b)  illustrates  positive  bias  (e.g.,  inventory  level 
after  starting  a  system  -fully  stocked).  Figure  1(c)  is  a  transient 
mean  process  which  damps  out.  Processes  (a)  ,  (b)  ,  and  (c)  each  appear 

to  be  approaching  "steady  state";  this  is  indicative  of 
initialization  bias  dying  out  as  run-length  increases.  The  process  in 
Figure  1(d)  has  not  yet  approached  steady  state  (and,  in  -fact,  may 
never ) . 

It  is  unrealistic  to  expect  that  4^  **  4  -for  all  i.  Hence,  we 
will  only  be  interested  in  detecting  significant  initialization  bias. 

The  tests  to  be  described  in  this  paper  can  be  motivated  in  the 
-following  ANOVA  sense:  We  partition  the  process  X,,...,Xn  into  two 
contiguous,  nan-overlapping  portions.  For  a  particular  realization  of 
the  process,  a  variance  estimate  based  solely  on  the  first  portion  of 
the  output  (or,  alternatively,  on  the  entirety  of  the  output)  is 
calculated.  A  variance  estimate  from  the  latter  portion  only  of  the 
output  is  also  calculated.  If  the  two  variance  estimates  are  deemed 
to  be  significantly  different,  then  we  reject  H^s  4^  *  4  for  all  i. 
CFor  the  transient  mean  processes  illustrated  in  Figures  1(a),  (b) , 
and  (c)  ,  we  would  expect  a  variance  estimate  formed  from  the  first 
portion  of  the  output  to  be  greater  than  an  estimate  from  the  latter 
portion.  3 

2.2.2  A  Test  Based  on  the  Maximum  Estimator 

With  the  intent  of  applying  the  theory  of  standardized  time 

series,  Schruben  (1982)  assumes  that  the  Y.  's  satisfy  the  mild 

1 

requirements  alluded  to  in  Section  2.1.  The  entire  -CX^  >  process  is 
then  standardized  into  one  <Tn(t)>  process. 
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Under  the  "stationarity"  hypothesis  HQt  u,.  ■  u  -for  all  1,  recall 

that  the  maximum  estimator  lor  variance  (when  the  number  of  batches 

b  ■  1)  is  given  by  V_  ,  s:  (c2/3)X2(3),  where  c e  lim  nVar(£X./n)  *= 

3,1  n-4®  i 

lim  nVar(£Y./n>  and  the  notation  "cM  is  read  “is  approximately 

n-4®  l 

distributed  as". 

Schruben  also  gives  a  variance  estimator  o  based  exclusively  on 

~2 

the  latter  portion  of  the  stochastic  process.  c  arises  from 

Fishman's  autoregressi ve  time  series  modelling  method;  Fishman 

2  2  ~2 

supposes  that  vc  c  x  (v)  ,  where  v  must  be  estimated.  Note  that  c 

and  V_  .  are  not  necessarily  independent. 

0,1 

~2 

If  we  nevertheless  assume  that  V_  and  e  are  independent,  then 

0,1 

F  s  V_  ,/er2  s?  F(3,v>,  (2-1) 

o  f  1 

the  F  distribution  with  3  and  v  degrees  of  freedom.  Let  f  be  a 
realization  of  the  random  variable  F.  If  we  assume  that  the  form  of 
the  initialization  bias  is  arbitrary  (negative)  Cpositi vel ,  then  we 
reject  HQ  at  the  a  level  if  f  <  f3fV>a/2  or  >  *3<Vfl_a/2 

<f  >  f3,v,l-«>  >  f3,v,l-a3’  Where  *3,v,V  iS  the  Upper  V  quantile 

of  the  F(3,v)  distribution  and  f  is  the  realization  of  F  resulting 

from  the  process  — X ^ ,-X^, • • • ,-X  • 

The  above  testing  procedure  appears  to  work  well  for  the  battery 

of  simulated  systems  studied  in  Schruben  (19B2).  However,  the  test 

may  not  perform  adequately  if  the  simulation  run  is  very  short  (since 

the  test  is  asymptotic)  or  if  initialization  bias  pervades  the  entire 

process  (in  which  case  the  experimenter  should  be  able  to  detect  the 

bias  visually).  There  are  also  a  number  of  problems  inherent  with  a  , 

^2 

including  the  fact  that  o  and  V_  .  are  not  necessarily  independent. 

O  g  1 

Schruben  gives  an  alternative  estimator  to  c  :  It  is  suggested 


.  AJ  ...1  -  A  A.*  .SI  » 
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that  the  proceiss  be  divided  into  b  »  2  batches,  each  of  size  m  *  n/2 
so  that  (assuming  the  batches  are  "approximately  independent"). 


h2(m-h2> 

.  \  s 

Kj (m-hj  >  S2 


3K  ^ (m-K  > 
-  __ 

mS2 

3K2(m-K2> 


F (3,3)  , 


(2-2) 


where  the  subscripts  1  and  2  refer  to  the  estimators  from  the  first 
and  second  batches,  respectively. 

Me  will  investigate  the  natural  generalization  of  this  test 
statistic  in  Section  3. 


2.2.3  A  Test  Based  on  the  Area  Estimator 

Schruben,  Singh,  and  Tierney  (1983)  works  with  a  weighted  form  of 

the  area  estimator  for  variance  to  test  HQ  vs.  H ^ :  u.  =  u(l-a.)  for 

all  i,  for  some  arbitrary,  pre— speci f i ed  constants  a^  ,  i*l,...,n. 

CThis  test  is  unrealistic  when  p  =  0;  in  this  case,  an  alternative  of 

the  form  H. :  u.  =  u  +  a.  could  be  used.}  The  authors  standardize  the 
1  l  l 

entire  output  series  and  find  that  (under  certain  strong  assumptions) 
the  most  powerful  test  is  to  reject  when  the  statistic 

2  3  skSk 

is  large,  where  c  =  a  —  a  .  They  give  several  examples  and 

K  k  k+1 

arguments  which  show  that  the  experimenter  can  offer  "reasonable" 

choices  for  the  c.  's. 

k 


Since  T  (t)  -*  as  n  -*  ® 
n  t  ’ 


2  =  I?=l  ckkSk  ”  *  £k=l  ckTn<k/n>  5  V"  »  ckBk/n‘ 


(2-3) 


B 


This  implies  that  I 


Z  :  Nor (0,ng  v) 


where 


(2-4) 


v  i  yn  ,  yn  .  c.c  [min(i/n,  j/n)  -  ij/n-"!. 

t>el  t<icl  1J 

Given  the  's  it  is  simple  to  explicitly  calculate  v  [See  Goldsman 
(1984)  for  additional  discussion  concerning  weighted  area  variance 
eet i mators. 3 

As  in  Section  2.2.2,  Schruben,  Singh,  and  Tierney  calculate 

^2  ^222 

another  variance  estimator  a  and  suppose  that  vc  z  a  \  (v).  The 

same  comments  and  caveats  as  before  still  apply  to  v  and  a  .  Further 

•'2 

supposing  that  Z  and  a  are  independent  yields: 


Z 

/  2  ,1/2 

(nc  v) 


Z/(ng2 

(c2/g 


v) 


1/2 


2}  1/2 


t  (v>  . 


Depending  on  whether  we  wish  to  test  for  arbitrary,  negative,  or 
positive  initialization  bias,  the  appropriate  test  should  be  performed 
as  in  Section  2.2.2.  The  Schruben,  Singh,  and  Tierney  test  procedure 
appears  to  work  well  for  the  examples  given  in  their  paper.  We  will 
generalize  this  procedure  in  the  next  section. 


3.  A  NEW  CLASS  OF  TESTS  FDR  INITIALIZATION  BIAS 

In  the  ensuing  discussion,  the  various  variance  estimators  given 

in  the  previous  section  will  be  used  to  construct  new  tests  for 

initialization  bias.  Divide  X.,...,X  into  b  adjacent  batches,  each 

1  n 

of  size  m.  Variance  estimators  based  on  the  first  b'  batches  will  be 
compared  to  the  corresponding  estimators  from  the  remaining  b-b ' 
batches.  This  comparison  is  to  be  accomplished  via  an  F  test. 

Using  the  notation  and  results  from  Section  2,  we  have  (under 


H^i  u  «■  u  <or  all  i  >  s 


0,b  • 


e  x^b'-l),  1  v  b"  \  b-1  (classical ) , 


2  2 

Oj  b'  -ex  (b  ’  )  ,  1  <  b  '  b-1  (a re?)  , 


^2^b’  :  c  %  <2b'-l),  1  <  b'  <  b— 1  (combined  classical-area) , 


2  2 

0-r  .  .  :  c  x  (3b),  1  <.  b  ’  b-1  (maximum). 


^4^-  ~ox  (4b  -1),  1  <  b'  <  b-1  (combined  cl  assi  cal  -max  i  mum) 


By  similar  reasoning. 


Q0,b-b-  s  «  Ii=b-*lC\,,n  -  bV  8-b-l  xj  ,J2  ••  »2*2<»-b->>. 


1 ,b-b * 


°l,b  ~  °lfb* 


2  2 

cr  x  (b-b  >. 


2,b-b '  ~  “o,b-b'  “l,b-b 


5-  a2X2  (2b-2b '-1)  , 


3,b-b ’ 


4,b-b 


Q-f  .  -  Q 

3  ,  b  3 ,  b 


2  2 

cr  X  (3b-3b'>, 


Q0,b-b'  +  °3,b-b’ 


'  tr'^x'*"  (4b— 4b '  — 1 ) 


To  condense  notation  a  bit,  define: 


p-1  for  k=0 
p  for  k=l 
2p-l  for  k=2 
3p  for  k=3 
4p-l  for  k=4 


Then  for  all  k 


J.  .  .  :  >  and  Q,  .  .  .r  c  %  (d  ..  )  . 

l:  ,  b  k,b  k,b-b  k,b-b 


Further,  define  tor  all  k  and  b: 


D  D 

, ,  k ,b '  .  .  *  k ,b-b 

k,b  d.  .  k,b-b  d.  .  .  . 

k,b  k,b-b 


Clearly,  the  V's  are  variance  estimators  calculated  from  the  first  b' 

-M 

batches,  and  the  V  's  are  the  analogous  estimators  from  the  remaining 


b— b '  batches.  Under  the  assumption  of  independent  batches,  Q  is 

K  ^  f  D 


i ndependent  of 


Q.  K  .  ,  »  for  all  k  ,k, 

k2,b-b  1 


We  could  thus  consider  test  statistics  of  the  form  (under  H^) : 


k^b 


k2,b-b 


F  (d.  U1,  d,  .  .  , )  • 

klfb  k2,b-b 


(3-1 ) 


One  could  attempt  to  use  a  test  statistic  of  the  form  V  .,/<r  »  but 

K  ^  9  D 

since  the  numerator  is  not  necessarily  independent  of  the  denominator, 
this  statistic  might  not  be  distributed  as  F(d  ,v) . 

Kj*D 

The  test  statistics  from  Section  2.2.2  are  easily  seen  to  be 
special  cases  of  the  above  CN.B.  (2—1)  and  (2-2)  have  m  =  n  and 
m  =  n/2,  respectively.! 

For  simplicity,  we  will  only  work  with  test  statistics  of  the 
form  (under  H^> : 


k  ,b  ' 


k  ,b-b ' 


F (d.  .  ,  ,  d  .  . ) 

k,b  k,b-b 


The  goal  now  is  to  find  that  combination  of  k,  b',  and  b  which,  in 


some  sense,  yields  the  "optimal"  test  statistic. 


POWER  CAL CUL AT  1  ON*.}  TOR  THE  NEW  TESTS 

ft  rwsonable  criterion  ■for  comparison  amonq  tests  (with  fixed 

✓ 

level  of  significance)  is  power.  Consider 

j  B  ECX  3  ■  u,  where  X  .  is  the  j-th  observation 
from  batch  i;  i«l,...,b  and  j«l,...,m. 

vs. 

Hts  i  “  4>»  where  the  a.  .'s  are  pre-specif ied. 

4  4  t  J  1  i  J  MJ 

We  assume  that  the  form  of  the  initialization  bias  under  H1  is 

negative.  (The  cases  of  arbitrary  and  positive  bias  are  similar.) 
Then  (3-1)  implies  that  we  must  reject  at  level  a  if 

^k  ,b  '  ^k  ,b- b  '  ^  ^d  ,,d  t,l— a* 

k , b  k |b~b 

We  give  analytic  results  for  the  cases  k  *=  0,  1,  and  2 
(classical,  area,  and  classical— area) .  Limited  Monte  Carlo  results 
for  k  =  3  and  4  will  be  given  in  Section  4.6. 


4,1  Classical  Batched  Means  Tests 

Suppose  that  ECX.  3  =  p(l-a.  .)  for  all  i,j.  Define: 

1  »  J  i  *  J 

Y  =  X.  .  +  Ua.  .  (so  ECY.  .3  =  u>  , 
i  »  J  i|J  i  ,  j  i,j 


.  =  —  Y 

i  ,  J  j  6p=l  i 


5  1  rJ 

i  ,  j  j  6p=i  ai  f, 


_  rb  — 

b  Li = l  ai,m’ 


c  =  lim  mVar(Y.  ) 
m->®  l  ,m 


Assuming  that  the  Y  's  satisfy  the  mild  requirements  alluded  to  in 

4  1  J 

Section  2.1,  a  direct  consequence  of  Theorem  21.1  of  Billingsley 
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(1V68) 


is  that  Y.  =  Nor  ( u ,  0Z/m> .  Hence,  Si,  e  ^  5"  T  ,  X 

i  ,m  i ,m  m  L J«1  i  ,  j 


Nor(p<l-a.  ),  /m)  .  This  leads  to  the  followinq: 

1  ,  m 


Theorem  i 


m 


Vb  ,<X  -X  >2  =  .V[b-1,  (uZm/.Z,T“  ,<5,  -a  >Z) 

t«i  =  l  i  ,m  n  *-  ‘•1*1  i  ,m  o 


2.  rb 


2  2 

where  x  (d,^)  is  the  noncentral  X  distribution  with  d  degrees  of 
freedom  and  noncentrality  parameter  £>. 

Proof i  See  the  Appendix.  ^ 

The  theorem  immediately  implies  thats 
°o,b-  5  ‘o,b->  and  °o,b-b-  *  'o.b-b-’’ 


where  l 


C»,b  ' 


.  <u2«,/.2>  yb',p  -  i~  P.  I,  J2  and 
£•1  =  1  L  i,m  b  £<j  =  l  j,m-' 


0,b-b 


s  (u2m/ff2)  yb  .  ,  ,  Qa.  -  .  1  ■ .  a  i 

£<i=b  +1  i  ,m  b-b  £*j=b  +1  j,m-' 


Hence,  Rob.  a  *  F(b  -1  ,b-b  -1,  »0>b  •.  «0,b-b  ■’  ’  "here 

F (d  , d_ , y  , y  )  is  the  doubly  noncentral  F  distribution  with  d  and  d 
degrees  of  freedom  and  respective  noncentrality  parameters  ^  and  y ^ 
The  power  of  the  test  is  therefore: 

PrfReject  HQ|H1  true!  =  PrtR^.  >  V-l  ,b-b  -1  ,  l-a*  * 

4.2  Area  Tests 

The  standardized  time  series  from  the  i— th  batch,  denoted  by 

<T.  (t)},  is  given  by: 

1  ,  m 


I mt 1 S . 

T.  _<t>  =  - V  t  1  ,  t  e  tO, ID. 


1  ,  m 


cVm 


Thus, 
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3  (J/m) 

i  ,m 


JS  .  j(X.  — X .  ) 

l ,  i  m  x  ,m  x O 

oVm  v  Vfti 


M  i  1 

T;  (j/m)  +  — y--  , 

i,m  crVm  * 


where  T  (j/m)  m  j(Y.  -Y .  .)/c^tn  and  c.  .  m  -j  (a  -a.  .),  with 

i,m  x  ,m  x,j  x,J  i  ,m  i  ,j 


the  Y  .  ’ s  and  a 
*  ,  J 


and  a.  's  defined  as  in  Section  4.1.  Hence, 

1  » J 

A.  =  ym  ,  jS  .  -  Vm  o  y1"  T.'  (j/m)  +  4  ,  c .  .. 

1  t-j“l  i,j  L j«l  x  ,m  t.j«l  1  , j 


Note  that  m<->  converges  to  a  Brownian  bridge  asymptotically  in  m. 
It  is  now  easy  to  see  Ccf .  Schruben  (1983)3  that  for  large  m: 


*  _r  -  *  ,  .  Cm 

lie  e  EC  A.  3  s-  u  >  .  .  c.  and 
1  1  t,j  =  1  1  ,  j 

<r2v*  =  Var(A.)  =  Var  fyin  c  ?m  ,  T.'  (j/m)}  5 

1  x  L  <-j  =  l  x  ,m  J 


m^-m  2 

IT  '  • 


So  as  m  becomes  large,  these  results  imply  that 


A  •  ?  #  a  a  a 

Ai  s  Nor(uc. ,ff“v.)  for  all  i.  Since  v  =  Vj  *  ...  =  vfa,  we  have: 


1  rb  £2  2,2r.  .  2.  2  *  rb  ,  *.2-» 

*  ii=l  Ai  *  *  *  Lbi  <M  /«r  v  )  2Li  =  1Cc.)  J  . 


2,  2*  rb  ,  *v2i 


Then 


Q  c*  #r 

l,b' 


■x^[b',  (u2/ff2v*)  J^jCc*)2]  and 


Ql,b-b'  '  ^2[b“b',  (u2/r2v*)  ^b.  +  1(c*)2). 


l,b' 


1 ,b-b ' 


*  F (b ' ,  b-b',  h.  .  e.  .  ..), 

1  fD  1 , D~D 


where  Ej  and  6^  ,  are  the  obvious  noncentral i ty  parameters. 


power  is  Pr<Rlb.  >  f b . ,b_b . , 


4.3  Classical-Area  Tests 


Cl  early. 


Q„  .  ,  =  Gk  .  ,  +  Q,  .  .  <r2x2(2b'-l,  E_  ,_  . )  , 

41  f  o  v  f  b  1 1  b  2  f  b 


where 
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'2,b- 


m  & 


0,b 


lfb 


Bi  mi  1  ar  1  y , 

°2,b-b '  "  Q0,b-b'  *  Ql,b-b' 

where 

*2,b-b '  3  *0,b-b '  +  *  1 ,b-b 


c2X2(2<b-b 


>-l 


*2,b-b') * 


So 


r  B  — JS  F(2b  -1 ,  2(b-b')-l, 

2  f  b  2?D“D 

2,b-b  ' 

which  has  power  PrCR^.  >  <2b • -! , 2 (b_b • > _* , !_a> • 


4.4  Weighted  Area  Tests 

As  an  aside,  we  describe  some  additional  tests  for  initialization 
bias  which  are  similar  to  those  in  Schruben,  Singh,  and  Tierney 
(1983).  Once  again,  test  H^s  wu  «=  u  for  all  i  vs.  HjS  =  u(l-a^) 

for  all  i.  We  construct  a  test  statistic  based  on  (2-3)  from  the 
first  b'  batches  of  the  Xj,...,Xn  process.  This  statistic  is  compared 

•» 

to  V.  ,  ,  in  order  to  test  for  bias. 

K  v  D“D 

By  (2-3)  and  (2-4) , 


3  Ii=l  Cii<Xn  '  V  *  N°r(0,nb.a2;b.)  , 


where  n,  ,  =  b'm  (b '  batches  of  m  X. 's  each),  c . =  a.  -  a.^«,  and 
b  i  lii+l’ 

2 


'  3  Ii  =  l  Ij  =  l  «=icj[min(i/nb.,j/nb.)  -  ij/nbJ  . 


Assuming  independent  batches,  Z  and  V 


k  ,b-b 


.  are  independent;  then 


under  HQ, 


~  2'  1/2 

7  /(  *W  w*  )1/2  -  Zb’/<nb’g  vb#)  . 

b'/<nb'Vb'Vk,b-b'>  '  (V*fb_b./e2)l/2  ‘  k-b“b'  ' 

Power  analysis  similar  to  that  from  subsections  4.  1  -  4.3  can  now  be 
performed. 


4.5  Analytical  Compa nson  of  lest* 

Mr  can  compare  the  power  o4  the  classical,  area,  and 
classical-area  test  •for  initialization  bias.  For  given  b,  b',  m,  and 
k  (■  0,  1,  or  2),  consider  the  alternative  hypothesis  H^s  EEX^  j -1  " 

u(l-a.  .),  i  *=  l,...,bt  j  =  By  previous  work,  the  statistic 

1  » 1 

R.  .  .  s  F (d.  ^.,d.  .  .  ,  ,  £.  c .  ...)  and  has  power 

k,b  k,b  T  k,b-b  k,b  k,b-b 

PrfR  .  >  f  .  ,  ,  3 . 

k’b  dk ,b ' ,dk ,b-b  * ’ 1  ° 

Since  tables  for  (singly  and  doubly)  nonventral  F  distributions 
are  not  readily  available,  we  approximate  the  distribution  of  R^  ^ ,  by 

using  the  familiar  result  Ccf.  Johnson  and  Kotz  (1970),  pg.  1973  that: 


F(V),  v2,  Vj,  V2)  '  cF(vlt  v2)  , 


where 


(v  +V  )v 
1  *1  2 


<W 


<V2+V2,V1  ’  1  wi+2vi 


,  and  v  = 


’W 

v2-2v2 


This  gives: 


k,b'  •'  (d, 


k,b  ^&k,b  )dk,b-b-  F[(dk,b<','&k,b,>  <dk,b-b  *ck,b-b') 
,b-b+ck,b-b)dk,b'  (_<dk,b  '+2^k,b  ' J  ’  (dk  ,b-b  -*2kk  ,b-b  '  *  a  ‘ 


So  the  power  is  approximately: 


f  F(vt ,v_>  >  -  f  .  ]  , 

l  12  c  dk,b ’ fdk,b-b • **  a  J 


where  c,  Vj,  and  v2  are  the  appropriate  quantities. 

Since  c,  v^,  and  «2  are  functions  of  ^ .  and  «k  . ,  they  are 
2  2 

functions  of  u  /v  ,  which  is  unknown.  If  we  are  willing  to  estimate 
2  2 

the  value  of  u  /o  ,  that  combination  of  b,  b',  m,  and  k  can  be  found 
which  maximizes  (*). 


ft  simple  example 

Suppose  that  m  =  100,  b  =  20,  and  H,:  ECX.  .3  =  u(i-a.  ),  where 
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a  ■  fl  -  — — -for  all  i,j.  l-igurec  2  through  7  give  plot*  of 
1  ,  j  *•  bm 

o  2 

<«*)  as  a  function  of  )i~ / v  for  b'  *  2 , 5 , B ,  1 0 ,  1 5 ,  and  IB, 

respectively.  Each  of  the  six  figures  illustrate  plots  for  the 

k  0,  1,  and  2  cases  with  the  level  a  *  0.05. 

Suppose  we  fix  b'  «=  2.  From  Figure  2,  it  is  easily  seen  that  for 
2  2 

all  values  of  u  /o  c  <0,51,  the  power  for  the  specific  alternative 
hypothesis  is  maximized  when  k  =  1  (area  test). 

A1 ternati vely,  fix  b'  =  B  and  consult  Figure  4.  We  see  that  the 

2  2 

power  (*)  is  maximized  for  u  /e  e  (0,.5)  by  k  *  O  (classical  test) 

2 

and  for  u  / o“~  «  [.5,51  by  k  =  2  (combined  cl assi cal  — area  test). 
Finally,  Figure  7  reveals  that  when  b'  «  IB,  (*)  is  maximized  for 
U2/*2  c  (0,51  by  k  =  2. 

The  above  results  are  to  be  expected  if  we  note  that  the 
classical  and  combined  cl assi cal -max i mum  variance  estimators  are  more 
sensitive  to  between  batch  changes  in  mean  than  is  the  area  variance 
estimator . 

4.6  Empirical  Compari  son  of  Tests;  An  Example 

We  empirically  compare  the  power  of  the  classical,  area, 

cl assi cal —area ,  maximum,  and  cl assi cal -max i mum  initialization  bias 

tests  for  the  following  example: 

Consider  the  stationary  first  order  autoregressi ve  process 

X.,.«.,X  ,  wi th 
1  n 

2 

X.  =  BX.  .  +  e.,  where  e.  ~  iid  Nor(0,l-fJ  ). 
i  l-li  i 

Let  Y  .  =  X  .  ♦  p(l-a.  .),  i=l,...,b;  J*l, - ,m.  We  test 

1  ,  J  »  » J  1  i  J 

l-L:  ELY.  .  1  =  u  for  all  i,j  against  H  :  ECY.  .3  =  u(l-a.  .).  Suppose 
O  1 , J  1  1 , J 


Lot  the 


I  / 

we  take  u  *■  1  and  m,  b,  and  a.  tk  in  the  previous  example. 

*  »  -) 

level  a  •  0.05. 

For  each  of  various  values  of  b'  and  AR  <  1 )  coefficient  B,  we  ran 
500  independent  Monte  Carlo  experiments;  we  applied  the  -five 
initialization  bias  tests  <k  »  0,1, 2, 3, 4)  to  each  experiment.  For 
given  b',  B,  and  k,  the  estimated  power  ■«  (number  of  experiments  (or 
which  Hq  is  r ejected ) /500.  Table  1  summarizes  these  results. 

Remarks: 

<1)  For  this  example,  the  maximum  test  (k  =  3)  or  the 
classical -maxi mum  test  <k  =  4)  appear  to  be  the  most  powerful. 

(2)  For  this  example,  it  turns  out  Ccf.  the  proof  of  Result  5-5  in 
Goldsman  (19B4)3  that  <r2  =  <l+B>/<i-B>;  sd  u2/<t2  *=  <  1-B  >  /  <  1  +  B>  .  We 
then  see  that  the  k  =  0,  1,  and  2  entries  in  Table  1  closely  match  the 
results  in  Figures  2  through  7. 

5.  CONCLUSIONS 

We  have  constructed  a  general  family  of  tests  for  detecting 
initialization  bias  in  a  simulated  process:  The  process  is  divided 
into  two  adjacent,  non— overl apping  portions.  A  variance  estimate  is 
then  calculated  from  each  of  the  two  portions.  If  these  estimates  are 
deemed  to  be  significantly  different,  then  bias  is  said  to  be 
present.  Of  course,  a  good  deal  of  the  simulation  literature  deals 
with  the  question  of  variance  estimation;  our  estimators  are 
primarily  rooted  in  the  recent  standardized  time  series  work,  but  use 
of  other  variance  estimators  is  strai ghtf orward . 

A  criterion  of  desirability  for  a  particular  initialization  bias 


IB 

test  ie  that  it  be  powerful .  ( Wp  were  able  to  derive  analytic  power 

results  for  the  classical,  area,  and  combined  cl assi cal —area  cases.) 
However,  as  the  examples  of  the  previous  section  show,  there  is  not 
necessarily  a  choice  of  k,  b,  and  b’  which  yields  the  most  powerful 
test  in  all  situations.  Therefore,  the  authors  are  currently 
addressing  this  problem  of  determining  the  best  k,  b,  and  b'  for 
general  processes. 
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Appendix 

We  prove  the  theorem  in  Section  4.1;  viz.,  if 


X,  -  Nor ( p(l-a .  ),  — ) , 

1  fni  i  ,m  tn 


,  ,  2  b 

,2  _  2  2r.  ,  pm  , 


i=»l  ’  0  i-1  * 


Proof ;  Consider  independent  U.  ~  Nor  (v^,!  ),  i«l,...,b.  There  exists  a 

nonrandom  orthogonal  bxb  matrix  H  whose  last  row  is 
_!/•>  -1/2  -1/2 

(b  “".b  ,...,b  ).  (An  example  of  such  a  matrix  is  easy  to 

T  _  _ 

construct.)  Define  Z  =  (Z Z.)  by  Z  =  UH  .  Then  Z,  **  /b  U.  ,  where 

l  D  ~  ~  b  b 

\  5  i  Ih  "i- 

Since  H  is  orthogonal,  ZZT  «  UHT(UHT)T  «  UHTHUT  »  UUT;  so 

rb  2  vb  2 
^i=lZi  "  ^i-1  Ui*  Thus* 


s2  s  l  (urub)2  -  l  U2  -  bU2  *  l  Z2  -  z2  -  h~l  z\. 

i-1  i-1  i-1  i-1 


T  T 

Since  H  is  a  nonrandom  matrix,  E[Z]  =  E[UH  ]  -  vH  .  The 
variance-covariance  matrix  of  Z  is  given  by: 


DU1  -  eK2~vH  )  (Z-vH')J 


T  T  T  T  T 

-  El(un  -vh  )  (yii  -vh  )] 

-  E[H(y-v)T(u-v)HT] 

-  HE[(y-v)T(U-v)]HT  -  HIt2HT, 


because  the  U^'s  are  Independent  Nor(v^,T  )  random  variables.  So 
D[Z]  -  It2. 

Since  each  is  a  linear  combination  of  normal  random  variables,  we 


have: 


Z  ~  Nor.  (vHT,It2). 

^  D  ~ 

I.e.,  the  Z^'s  are  independent  normal  random  variables  each  with  variance 
2 

T  .  Thus , 


2 

T 


b-1 


I 

i=l 


x2<b-l,6). 


2 

where  6  is  computed  as  follows:  It  Q  ~  x  then  E [ Q ]  =  b-1+6  [see 

Rohatgi  (1976),  pg.  315].  This  yields: 


ElS‘/f  ) 


-  -J  {ei  }  r;i  -  bEiubn 


{  ]  IVarCUj)  +  (EU^2]  -  biVar(Ub>  +  (El^)2]} 
1  i-1 


— 2  {bi2  +  \ 

T  i-1 


2-2  -  1 
-  T  -  bv,  },  where  v,  =  — 
b  b  b 


b 

l  \ 

i-1 


-  b-1 


b 

l 

i-1 


<VV 


2 


So 


.  1  r  (  -  .2 

5  =  "7  l  ^vi"vb^  ’ 
t  i-1  1 

The  proof  is  completed  if  we  identify  o  /m  with  t  ,  X.  with  U 

i,m 

and  u(l-a4  )  with  v. ,  i  =  l,...,b. 

l ,m  l  /  / 

An  alternative  proof  is  given  in  Appendix  A. 3  of  Goldsman  (1984). 


(d)  a  transient  mean  process 
which  does  not  appear 
to  be  approaching 
steady  state 


Figure 


V.  Various  transient  mean  processes  (drawn  as  continuous 
functions  for  clarity). 


k  •  1  (ares  test) 
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7  1.54  2.31  3.08 

RATIO 

The  power  (*)  as  a  function  of  the  ratio 
example  described  in  the  text  for  b'  *  5 
N.B.  The  ordinate  axis  runs  from  0.00  to 


- 1 

3.8 

and  k 

0.10. 


*  *  2  (carsined  classical-ires  test 


?  ? 

The  power  (*)  as  a  function  of  the  ravo  ■,  /c  for  the 
example  described  in  the  text  for  b'  *  8  and  k  *  0,1,2. 


.--V-V 
'  V.vIV 


Ficjre 


POWER 


a  =  -0.5 
k  =  0 
k  =  1 
k  =  2 
k  =  3 
k  =  4 

a  =  0.0 


=  c 

0 

o 

1  J 

i  u 

.000 

.000 

.396 

.970 

.992 

.270 

.058 

.060 

.092 

.088 

.066 

.066 

.000 

.000 

.626 

.996 

1.000 

.812 

.056 

.082 

.116 

.086 

.108 

.096 

.000 

.000 

.710 

.996 

1.000 

.948 

k  =  0  .000 

.002 

.260 

.658 

.762 

.170 

k  =  1  .070 

.076 

.050 

.064 

.044 

.064 

k  =  2  .000 

.005 

.304 

.732 

.886 

.422 

k  =  3  .090 

.108 

.076 

.098 

.110 

.080 

k  =  4  .000 

.012 

.368 

.742 

.942 

.652 

=  0.5 

k  =  0  .006 

.016 

.164 

.312 

.370 

.092 

k  =  1  .050 

.054 

.038 

.044 

.060 

.068 

k  =  2  .006 

.022 

.166 

.290 

.436 

.156 

k  =  3  .066 

.108 

.132 

.122 

.124 

.110 

k  =  4  .010 

.052 

.208 

.364 

.532 

.322 

Table  1 

Estimated  power  of  tests  for  initialization  bias  for  a  shifted  AR 
process.  For  given  b'  and  a,  the  five  table  entries  are  base 
on  the  same  500  (#  =  1000)  independent  experiments.  The  level 
of  significance  =  0.05.  See  the  text  for  details. 
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